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Abstract 

The problem of quasistatic and rate-independent evolution of elastic-plastic-brittle 
delamination at small strains is considered. Delamination processes for linear elastic 
bodies glued by an adhesive to each other or to a rigid outer surface are studied. 
The energy amounts dissipated in fracture Mode I (opening) and Mode II (shear) 
at an interface may be different. A concept of internal parameters is used here on 
the delaminating interfaces, involving a couple of scalar damage variable and a plastic 
tangential slip with kinematic-type hardening. The so-called energetic solution concept 
is employed. An inelastic process at an interface is devised in such a way that the 
dissipated energy depends only on the rates of internal parameters and therefore the 
model is associative. A fully implicit time discretization is combined with a spatial 
discretization of elastic bodies by the BEM to solve the delamination problem. The 
BEM is used in the solution of the respective boundary value problems, for each 
subdomain separately, to compute the corresponding total potential energy. Sample 
problems are analysed by a collocation BEM code to illustrate the capabilities of the 
numerical procedure developed. 
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1. Introduction 



Applications of layered structures are numerous and continuously increased, 
an example being the massive use of composite materials in aeronautical in- 
dustry at present. Usually the interfaces between these rather bulk laminates 
consist of very thin adhesive layers. For efficient computations, these adhesive 
layers may be approximated by zero thickness interface layers. There are many 
situations where an adhesive layer is found to be partially or fully damaged. 
This process is frequently referred to as delamination or debonding of adjacent 
material laminas. In this work the description of the damage is based on a scalar 



damage quantity (variable, cf. [IJ]), which is defined at interfaces and takes val- 



ues from the interval [0,1], with zero value meaning no adhesion due to the 
total damage of the adhesive while the unit value meaning complete operation 
of the adhesive without any damage. During a damage evolution the damage 
variable decays in time, and it is assumed that a specific amount of energy has 
to be released (dissipated). This simplified approach, motivated essentially by 



GriflSth jig , is often inadequate as it is observed experimentally that consider- 



ably more energy is usually needed to perform delamination in shear Mode II 
than in opening Mode I [TqI . 24 , i^j . Motivated by the microscopical idea of 
interface plasticity (23. |42|. an extra inelastic parameter is introduced 3^, 37 1, 



which describes some plastic slip that may occur in the tangent direction of an 
interface before its debonding. 

An alternative approach to model fracture-mode-scnsitivc delamination uses 
only the delamination variable but makes the dissipated energy directly depen- 
dent on the so-called fracture mode mixity angle, cf. (fTO l-lfTTj) below. This 
approach has frequently been used in engineering models 40|, |41| but, it does 
not seem amenable to a rigorous mathematical analysis. In the present work we 
consider the delamination as a unidirectional process, i.e. no healing (or recon- 
struction) of adhesive is allowed, which covers most of engineering applications. 

The goal of this article is to present and analyse from an engineering as well 
as numerical implementation viewpoint some basic features of the delamination 
model devised in [s^ [s^l with different dissipated energies in Modes I and II. 
In particular, in Section [2] we briefly present the energetic approach employed. 
In Section [3l we concisely introduce the present model, while some engineering 
insight on this model is provided. Then, in Section 21 the numerical imple- 
mentation of the model, is presented. Finally, in Section [5l two-dimensional 
simulations are developed, showing that the model is suitable for solving realis- 
tic problems of delamination between elastic layers. 



2. Theoretical background 

2.1. Problem definition 

Let us consider an assemblage of elastic bodies, each of them defined by 
a reference domain fli [i = 1, N), with the Lipschitz boundary Vi = 9f2i, see 
Fig. [TJ We denote by Tij = d^i n d^j the (possibly empty) interface boundary 
between f2i and Slj (i, j = 1,...,A), which may undergo delamination. We 
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Figure 1: Schematic illustration of the geometry and notation for a two-dimensional case of 
two bonded subdomains, i.e. N=2. 



also consider possible delamination on some parts of the outer boundary Foi, 
which is assumed to be in Signorini elastic contact with a fixed rigid surface, 
see Fig. [1] The union of these parts is denoted as Fg = Ui<i<Ar ^oi- We will 
denote Fc := {Ji<i<j<N Tij UFq. We assume that the rest of the outer boundary 
9f2 is the union of two disjoint subsets Fo and F^, where Dirichlet (prescribed 
displacements ud = uu{t)) and Neumann boundary conditions (prescribed trac- 
tions pn = PN{t)) are imposed, respectively. For the sake of simplicity of the 
following considerations, vanishing tractions = will be considered here- 
inafter, except for Sections 14.21 14.31 and 15.1.21 The intersection of the closures 
of Fc and F^ is assumed to be the empty set, i.e. Fc n Fd=0. Any Tij is con- 
sidered as an infinitely thin adhesive layer, represented by springs distributed 
continuously, similarly to the Winkler spring model, with distinct normal and 
tangential elastic stiffnesses of values ranging from zero to infinity. Both the 
elastic subdomains and the adhesive layers are assumed to store energy, which is 
given by a stored energy functional S'{t, u, z) a function of time i, the displace- 
ments u and the inelastic (internal) parameters collected in z. It is considered 
that two elastic subdomains VLi and f2j, may debond along the interface F^. 
During this process the material of the adhesive can be damaged and plastified. 
The onset and growth of the damage and plastification, represented by the z 
variables, does not depend on some internal time scale and therefore the process 
is considered as rate-independent. The damage and plastification of the adhesive 
layer are accompanied by a release of stored energy. The dissipation potential 
^(z), with z := 4|, for a rate-independent process can be represented by a 



degree- 1 homogeneous functional [31|. The processes described in this work are 
assumed to be quasistatic, i.e. no inertia effects are taken into account. The 
rate-independent evolution we have in mind is governed by the following sys- 
tem of doubly nonlinear degenerate abstract static/evolution inclusions., referred 
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sometimes as Biot's equations generalizing the original work [5|, |6[: 

du'g'it, u, z) 9 and + d^.S'it, u, z) 9 0, (1) 

where the symbol "9" refers to a (partial) subdifferential, relying on that ^(•), 
^(i, -jz), and S'{t,u,-) are convex functionals. The first optimality condition 
of Eq. ([1]) represents the minimum, energy principle^ while the latter one, the 
minimum dissipation potential principle i36| . 

For the sake of simplicity, throughout this work, we will restrict ourselves to 
the two-dimensional case, i.e. fi^ C will be planar domains, i ~ 1, ...,iV, and 
Fi, will be one-dimensional surfaces. 



2.2. Energetic solutions 

A fruitful concept of a certain weak solution to the doubly nonlinear inclusion 
with degree-1 homogeneous dissipation potential called energetic solutions. 



was developed by Mielke et al. [32|, |33[. In the convex case, this concept is 



essentially equivalent to conventional weak-solution concept, while in our case 
where ^{t, •, •) is non-convex this concept represents a certain generalization; cf. 



27[ for a survey on the concept of energetic solutions and [28[ for comparison 
with other concepts. 

The process {u{t),z{t)), t <S [0,r] is called an energetic solution to the 
initial- value problem ([T]), if it satisfies the following three conditions: 

(i) The energy equality: 

S[T,u{T),z{T)) + DissK(2;[0,T]) 



stored energy energy dissipated 

at time t = T during [0, T] 



;'(t,U,z)t + <f(0,Wo,Zo) , (2) 







work done by stored energy 

mechanical load time t = 



N 

where DissK(z; [0, T]) sup ^ 7^(z(^,) - z(<,_i)), (3) 



with the supremum taken over all partitions 
< to < <i < ••■ < tN-l <tN <T. 

(ii) Stability inequality for any t G [0,T]: 

S'{t,u,z) < S'{t,u,z) + ^{z-z) forany(w, z), (4) 

(iii) The initial conditions: u(0) = uq and z(0) = zg. 

In Eq. ([2]), tf/ is the partial derivative of S' with respect to time t. 
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3. Model of interface damage and plasticity 



In this section we present the specific model adopted, in order to simulate 
the nonlinear inelastic behaviour of an adhesive layer, by defining a suitable 
stored energy functional as well as a dissipation potential. The present plastic- 
type model with kinematic-type hardening [igI . Isoj for the delamination prob- 
lem, was devised essentially in (36| without any mathematical or computational 
justification, and further scrutinized in j37| . Beside the displacement u, two 
internal parameters are used in order to describe the nonlinear behaviour of the 
adhesive: the damage variable C and the plastic tangential slip variable tt, which 
together constitute the pair of inelastic variables z = {C,t^)- 

3.1. Stored energy 

Stored energy (S includes the clastic bulk contribution and the additional 
adhesive-surface contribution: 



S'{t,U,z) = S'cl{t,u) + (fadh(u, 2^) 



(5) 



with 



(?el(t,u) =< 



Ez=i/n. Qe(M):e(M)x, 



if u\to = uu{t), 
elsewhere. 



where C,; is the elastic moduli tensor in r2, , and 



(6) 



rc 



oo 



c(yH:+y(H,-)^ 

if < C < 1 and 
Mn > on Tc, 
elsewhere, 



(7) 



where Kn > and Kt > are the phenomenological elastic constants describing 
the stiffnesses of the linearly elastically responding adhesive in the normal and 
tangential directions, respectively, > is the so-called factor of influence 
of damage M = MnJ^ + Hfr with = and |Ju|t = M-t, ly 

and T being unit normal and tangential vectors to Fq, and ds is the tangential 
derivative defined on r^). For Fq the outward normal v is typically taken. 
Constant parameter stands for the plastic modulus of kinematic hardening. 
Here we used the notation |m] for the differences of displacements from both 
sides of Fc. We also assume r > 1. The last term in Eq. ([7]), although bearing a 
physical interpretation 0, is here introduced mainly for mathematical reasons 
in order to facilitate a proof of convergence; for further details see also [37l |. 
but in specific simulations one may expect reasonable numerical results even if 
this term is neglected by setting kq = 0. The constraint > in Eg. ([7]) is 
actually the Signorini non-penetration condition of unilateral contact [22|. 
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3.2. Dissipation potential 

A suitable functional for the dissipation potential, describing both inelastic 
processes of damage and plastic slip in the adhesive layer, and actually being a 
degree-1 homogeneous functional, is defined as: 

Cic |C| + 0't,yicld|7i'|S 

^(i)=^(C,^):=r^^ ifC<Oa.e.onre, 

otherwise. 

Parameter Gj^ > is the minimal energy required for complete damage (debond- 
ing) of a unit area of the interface. In particular, we assume it represents the 
interface fracture energy in Mode I. Parameter crt^yidd > is the interface yield 
shear stress for initiation of tangential plastic slip along the interface. The con- 
straint C < in ([8]) makes the evolution of C irreversible, i.e. the model does not 
permit healing, which means that a debond appeared at some point can not be 
restored. 

Note that, except trivial case when ud is constant in time, (f/ in ^ would 
not be well defined. One way how to avoid this drawback, well consistent with 
BEM, is to restrict the displacement only on Fc, assuming that Tc and To are 
not touching each other. The restricted displacement u|rc will be denoted by 
Uc', in fact, in the case of Fy. it is a couple of traces of u from both sides of 
Tij. As u does not occur in Eq. ([5]) and thus it is fully nondissipative, Eq. Q 
implies that u minimizes (o'{t, •, z) and thus, in fact, Uc and z determines m at a 
given time t. Thus, S" can be considered as a function of Uc instead of u, which 
makes <S'^{t,Uc, z) well defined if Ud is smooth in time. On the other hand, we 
will not distinguish between [[u] and [uc] . We will use this convention through 
the rest of this article. 

3.3. Engineering analysis of the traction-relative displacement law 

In the case of a linear elastic-brittle interface model (4^ |4l| , the interface 
failure criterion is connected to the energy release rate (ERR) concept. It can 
be shown [2^ that the energy stored in the adhesive at the crack tip equals 
the ERR of a mixed mode crack propagating along a linear elastic interface, and 
can be evaluated as: 

G = Gi + G„ = ^ + ^. (9) 

The so-called fracture mode mixity angles^ denoted as jpa, ipu, or ■0^, can be 
defined in terms of ERR as, 

tan2V'G = ^, (10) 
as well as in terms of relative displacements and tractions, respectively, 

tan t/Ju = IT- ir and tani/'cr = — = — riT- (H) 
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Thus, the fohowing relations hold: 

I tan-0(j| = . /— tan-0G and \ ta.n ipu] ~ \ — tan ipc- (12) 
V '^n V '^t 

It is assumed that a crack propagates if the ERR G reaches the fracture energy 
Gr, that means: 



Gc. (13) 



A strong dependence of Gc on the fracture mode mixity has been observed 
in extensive experiments [l , 12 , 1^, 24 1 . In accordance with other experimental 



observations l21|,|43, the associated plastic zones in the adjacent bulk, near the 
crack tip, are larger in Mode II than in Mode I and these plastic phenomena are 
localized in a relatively narrow plastic zone in the bulk in the interface vicinity. 
In order to provide a better representation of these experimental results, a 
plastic tangential slip variable tt has been introduced at the interface, which 
allows us, firstly, to distinguish between fracture Mode I and II in the sense 
that some additional dissipated energy is associated to interface fracture in 
Mode II, and secondly to simulate these narrow plastic zones. In such a case 
we can model an inelastic behaviour in the tangential response of the interface, 
while the response in the normal direction remains linear elastic, as shown in 
Figure [5] An engineering insight into the present interface constitutive law can 
be summarized by the two conditions which activate the two inelastic processes 



included in the formulation 371. The first one is the activation criterion for 



damage initiation which, for the case of ko=0, reads as 

l{^nlujl + .,{luj-n)')=G,^, (14) 

where the left hand side represents the elastic energy stored in the adhesive. 
The second one concerns the evolution of tt which is triggered when \at — Kj^Trj 
reaches the activation threshold Ct. yield, and then, 

|C'«t([w]j-7r) - Kjj7r| (Tt^yiold- (15) 

A more detailed analysis of the model may be found in [37| . The model produces 
the desired results if 

^V2KtG,, < at,yicid < V2KtG,,. (16) 

The upper bound of yield stress is necessary for making possible to initiate 
plastic slip before the total interface damage, while the lower one is required to 
avoid plastic slip evolution at some point which has already been debonded. 

Thus, the ERR of a mixed mode crack for the present model is defined by 
the: 

G = — 7: 1 7: 1- o-t,yicid TT H — , (17) 
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(b) Mt^yiold Mt,crit Mt 



Figure 2: Schematic illustration of the traction-relative displacement law in the model, (a) 
pure normal (opening) mode and (b) pure tangential (shear) mode, considering f o = 1 and 
ttq = 0. Contribution of the delamination-gradiont term is neglected, i.e. «;o=0. 
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where it may be seen that, referring to Eq. ERR is here augmented by terms 
concerning inelastic slip tt. 

In the following we will try to determine the dependence of Gc on the fracture 
mode mixity angle ■0„ similarly as in [i^. To accomplish this task, first we 
eliminate the plastic slip tt, which in our kinematic-hardening model may be 
written, for tt > 0, as: 

Substituting Eq. into the damage initiation criterion of Eq. , leads to 
the relation, 

\{'^'i"tM^Sm,^'-^f)-o,, (19) 

when some interface plasticity occurs, i.e. \u\t > ^^"^f^- 

In a similar way as in Eq. (|13p. which is valid if no plasticity has occurred, 
Eq. (|T9| defines the relation between the two components of the relative dis- 
placement at the crack tip leading to the crack growth, if some plasticity has 
already appeared. This relation can be written in a parameterized form through 
the use of a parametric angle (/>, as: 



2G, 



■ cost 



26*1^ Kt -I- Kh . , CTt yield .„„^ 

sm0 , (20) 



Cf, yield 



for arcsin— ^iS== < < f. Before plasticity occurs, i.e. for < < 
arcsin , the analogous parameterization writes as 



2G,. 



■ COS( 



Ht = V^^i^-^' (21) 



and angle cj) coincides with the fracture mode mixity angle ^pc defined in Eq. (jlOp . 
Parameterization of Eq. ((20|). defines an ellipse whose center is at the point 



(0, — ) , which continuously switches from the ellipse with the center at the 
origin of coordinates Eq. (PT|) , which corresponds to a state of zero plasticity. 

The relation Gc ~ Gc{ipu), for the case of non-zero interface plasticity, can 
be obtained by substitution of Eqs. ([201), ([19]) and (HH) into Eq. dTT]), leading 
after some algebra to: 

G,, f 1 + ^ sin^ A %M ^ Ge(0). (22) 
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Finally, finding the relation between the angles (j) and ipu 



. I I '^n '*t+'*H , , / '^n CT*, yield 1 

tanV'„ = W tan(/)-./— '-^ -, (23) 

we obtain the desired relation (j) = (l){tpu) to be substituted into Eq. p2|). How- 
ever, an explicit relation oiGdipu) is rather cumbersome. Nevertheless, accord- 
ing to plots presented in (isj the functional dependence of Gdipu) qualitatively 
represents the expected behaviour in view of the previous experimental results 



4. Numerical implementation 

The theoretical framework, briefly presented up to this point, provides an im- 
plementable and efficient numerical scheme. An emerging global minimization 
problem, inherent in Eq. (j4|) may be defined by an implicit time discretization. 
By discretizing the time incremental formulation in space by some appropriate 
method, the problem may be casted in a standard algebraic form. Since the 
problem may be (and here is) formulated on the boundary, as all nonlinear 
processes considered occur exclusively on the boundaries, Fc only, a boundary 
element method seems to be a natural approach especially if the bulk equations 
can efficiently be solved, which is, in particular, the case of isotropic linear elas- 
tic materials considered in this article. Such a formulation was developed in 



37[, using the collocation BEM but without providing a thorough description 
of the numerical implementation. A related symmetric Galerkin SGBEM for- 
mulation can be found in [43^ and a FEM implementation in [s^. Preliminary 
comparison with the SGBEM formulation has shown an excellent agreement in 
a few specific case studies. An advantage of the present approach with respect 
to a related FEM approach is that no bulk discretization is required here, and 
in the analysis and optimization procedures we directly work with a relatively 
small number of variables associated to boundaries in particular to Tq. 



Minimization problem 

Making an implicit time discretization by adopting, for simplicity, an equidis- 
tant partition of [0,T] with a fixed time-step r > 0, assuming T/r <E N, Eq. (|4]) 
leads to a recursive minimization problem: 

minimize ,^^{uc, z) ~ S{kT^ Uc, z) + [%{z—z^^^) 
subject to Biuc>0, 0<C<C''~^, 

to be solved successively for fc = 1, ...,T/t, starting from uq and zq. Operator 
Bj represents the non-penetration Signorini conditions, while the further const- 
raint in ()24p refers to non-negativity and irreversibility of damage parameter 
evolution. According to the convention of Section [3.2l only Uc, the displacement 
at interfaces (or contact zones), appears in Eq. p4p making clear that only this 
part of the displacement field is a minimizcr of the problem. 
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Table 1: Pseudocode of the alternate minimization algorithm 



(1) Set j = and C° = C''"^ 

(2) Repeat 

(a) Setj=j + 1 

(b) Solve for Uc and tt^ : 

minimize (uj,7r^) i-> + M{z^—z'^~^) 

subject to BjUc>0 

(c) Solve for C^ : 

minimize H> (g'{t, itj, z^) + ^{z^-z''-'^) 
subject to C*""^ > C-' > 

(d) If II (J-C^'-i II < e exit loop 

(3) Set = and z'' = z^ 



We denote by {u'^.^z'^) some (generally not unique) solution to the problem 
(El- 

In order to numerically solve the emerging minimization problems (|24p . we 
have utilized and test in this work several algorithms, such as the L-BFGS-B 0] 
for general large scale simply bounded problems, the GLPK routines for linear 
programming problems [25| as well as a conjugate gradient based algorithm with 
constraints, see Uj, for solving quadratic programming problems. Notably, 
the minimization problem appears to have an Li-type non-smooth term with 
respect to the plastic tangent slip variable tt. sec Eq. ([8]). In order to overcome 
this difficulty we take advantage of gradient projection algorithm presented in 
[l3| for such kind of non-smoothness. 

4-1.1. Alternate minimization algorithm 

The functional .^'^ in Eq. ([24| is not convex and as such leads to a difficult 
minimization problem. In order to overcome this difficulty we utilize a special 
technique, originally proposed in Q, called as alternate minimization algorithm 
(AMA). The AMA procedure, in our case, consists in splitting the original non- 
convex minimization problem to two distinct convex problems with respect to 
the kinematical variables (u, tt) and to damage variable Cj respectively. Con- 
vergence is succeeded through an iterative procedure by alternation of this two 
convex problems. A flowchart of AMA may be seen in Table [TJ It is worth 
mentioning that the individual sub-problems emerging by using such alterna- 
tion consist of a nonsmooth quadratic programming problem, step (2-b), and a 
linear programming problem, step (2-c) of Table [U respectively, for which we 
may use appropriate specialized algorithms such as those mentioned above. 

4-1.2. Back-tracking technique 

The above AMA procedure does not necessarily lead to a globally minimiz- 
ing solution which is, however, one of the main ingredient behind the energetic- 
solution concept, as shown in Section[2] In order to execute the global minimiza- 
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Table 2: Pseudocode of energy-based backtracking algorithm 



(1) Set k = 1 and C° = Co 

(2) Repeat 

(a) Determine C,^ using the alternating minimization 
algorithm for time tk and the initial value C° 

(b) Set C° = C'' 

(c) If the two-sided energy estimate holds: 

set k = k + 1 

(d) Else set A: = fc - 1 

(e) Until k > T/t 



tion more successfully at particular time levels we use heuristic back-tracking 
algorithm (BTA), devised and tested on such sort of problems in 0,14; 31, 36, 37 1. 



The BTA technique is based on checking a two-sided energy estimate, the inte- 
gral expression in Table [21 where also some pseudo-code of BTA is given. This 
two-sided inequality has been constructed by use of the energy stability condi- 
tion Eq. (|4]) and a full deduction of it can be found in (2^, [30[ • The upper 
and lower energy estimates are given as time integrals of the power while (mc, z) 
and {Uc,z) are piecewise constant interpolants in time defined by 

Uc{t) = for t e ((fc-l)T, fcr] , (25a) 

u^{t) = u^-i for t e [(fc-l)T, fcr) . (25b) 

Similar notation concerns also z and z. A thorough deduction of the boundary 
forms for these integrals of power, amenable into the boundary element context, 
is given in Sections 14.21 and [4.31 Although there is no proof that BTA converges 
to the global minimum, definitely it leads to solutions of lower energy than those 
obtained if we used AMA only. 



4-2. Boundary element method 

The boundary element method is closely related to the map between the 
prescribed boundary conditions in displacements or tractions and the unknown 
boundary displacements or tractions. In pure Dirichlet and Neumann boundary- 
value problems (BVPs), these maps are called Steklov-Poincare and Poincare- 
Steklov maps [lo, IS] , respectively, and BEM can be considered as an approach 
to discretize these maps. In the present computational procedure, the role of 
the BEM analysis, applied to each subdomain $7^ separately (which, in fact, 
makes this problem very suitable for parallel computers), is to solve the corre- 
sponding BVPs on each 51,^. For this goal, we numerically solve the Somigliana 
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displacement identity [35l |38| 



X 



(26) 



where ^ G = dfli and u'l^ix) and pl^ix) denote the m-component of the 
displacement and traction vector, respectively. The superscript i used in this 
section refers to domain Vli in difference to previous and next sections where 
it denotes time step. The weakly singular integral kernel J7^;(a;,^), two-point 
tensor field, given by the Kelvin fundamental solution (free-space Green's func- 
tion) represents the displacement at x in the m-direction originated by a unit 
point force at ^ in the /-direction in the unbounded elastic medium whose ma- 
terial properties coincide with those of D,i . The strongly singular integral kernel 
T^j(a;,^), two-point tensor field, represents the corresponding tractions at x in 
the m-direction. The coefficient-tensor c^/(C) of the free-term is a function of 
the local geometry of the boundary at ^, and may be evaluated by a closed 
analytical formula for isotropic elastic solids [2^ . The symbol ^ in Eq. (|26p 
stands for the Cauchy principal value of an integral. 

Consider a discretization of the boundary F^ by a boundary clement mesh, 
which is also used to define a suitable discretization of boundary displacements 
u^{x) and tractions p^{x) by interpolations of their nodal values. By impos- 
ing (collocating) the Somigliana identity (^5)) at all boundary nodes (called 
collocation points) we set the BEM system of linear equations for F^ . The so- 
lution of this system defines the unknown nodal values of displacements and 
tractions along F^ representing a part of arrays of all nodal values denoted as 
u* and p*, respectively. The arrays and p' also include the known nodal 
values of displacements and tractions along Fi given by the prescribed bound- 
ary conditions. The BEM system obtained from Eq. ([26]) is usually written as 
H*u* = G'p* [3^ . In our computer implementation of BEM, we employ straight 
elements with continuous and piecewise linear interpolation for displacements 
and possibly discontinuous piecewise linear interpolation for tractions. 

Then, to compute an approximation of the elastic energy, (foi from Eq. ([5]), 
stored in each bulk fi^ , by using the obtained approximations of boundary dis- 
placements ul^ and of the corresponding boundary tractions along F^, we 
utilize the following general relation [l^l, neglecting body forces: 




(27) 



while the corresponding total potential energy is 




(28) 
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Figure 3: Solution of a mixed BVP P for a single elastic domain given as a superposition of 
the solutions of the three sub-problems P'~^,P" and P^. 

Notice that both Eq. (j27p and (j28p provide pure boundary expressions of energy. 
If (foi in Eq. ([5]) is replaced by the sum of potential energies IIq. for all sub- 
domains, the functional S defined by Eq. ([5]) will represent the total potential 
energy of the whole problem. 

We will also need to compute integrals of time derivatives of energy, ap- 
pearing in Table [21 were it is presupposed that displacements on the contact 
boundary part Uc are defined in time steps k and fc— 1. In order to do such 
calculations we need to separate the problem into three different sub-problems, 
in each of them either the contact or prescribed Dirichlet or Neumann data 
are defined on the boundary. This separation to sub-problems may also serve 
to express and solve the minimization problem, considering the integral on Tq 
only. 

^.3. Boundary forms of the total potential energy for a single domain 

Consider the BVP for a sub-domain fi^. In this section we will omit index i, 
for the sake of simplicity. Let and p,,, respectively, denote the displacement 
and traction solutions of this BVP restricted to F^, 77 = C, D and N, e.g. 
Uc = u\tc and p^ = p\ru- We assume here a mixed-type operator ^ which 
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formally assigns (pc,Pd,Wn) to the known boundary data (mc,Ud,Pn) of the 
original BVP P shown in Fig. [3l and may be expressed using the following 
block structure as: 



(29) 



The columns of the aforementioned block operator ^ are associated to the sub- 
problems defined in Fig. |31 The displacement solution of a subproblem P^ 
is denoted as w'' . From the principle of superposition the displacement solution 
of P may be reconstructed by the sum: 

u = u'' +u° + , (30) 

The total potential energy for the mixed type BVP P can be written in an 
expanded form as, 

n(t, w)=-;r /ucPcS+i [ UdPdS-^ I u^p^S. (31) 
^Jrc ^Jto ^JTt, 

By substituting the unknown data for the problem P from Eq. ()29p the total 
potential energy writes as 



n(^,.c) = i(/^ 



Uc^cqUc 9 
Tc 

Fd Fd Fd 



V^Ji^cUc / Pn^ND^d S- / Pn^nnPn S 1 . (32) 
Fn -^Fn -^Fn / 

From Eg. ((32)) it is clear, that since Ux^it) and PN(i) are known, the total 
potential energy is in fact a function of the contact displacement Uq only, in 
addition to be a function of time t. We further modify Eq. ([5^ . in order to 
hold the unknown variables on Fc only, by utilizing the second Betti reciprocity 
relation between the elastic solutions of P"^ and P^ , 



Pn.^nc"c S = - / Uc./#cnPn S (33) 
Fn Fc 

as well as between the solutions of P'^ and P°, 

/ Uo^y^'ocUc S = / Uc.'^cdWd S. (34) 

JFd "'Fc 



15 



Then, by substituting Eqs. 1^ and ^ into Eq. ([5^ . 



Tc 

+ - /md-^ddWdS + - Uo^dnPnS 



-i /pn^-^ndWd S-i /pn.^nnPnS. (35) 

The partial time derivative of the total potential energy expression in Eq. psp 
writes as 



McI-^cdWd+^cnPn) S 
1 /■ ^ „ 1 



-- Pn^noUdS-^ /pn^nnPnS (36) 



where the bar with dot denotes the time derivative of the expression below the 
bar. The integral corresponding to that on the left-hand side in the two-sided 
inequality in Table [2J can be evaluated using Eq. p6p as, 



J(k-1)T JTc 

- I W^(^cdUo"^+.^cnPn"^) S 
JTc 



/ pI-^^^pI S+i /" p^~'^nnPn~' S, (37) 
and similarly for the integral on the right-hand side of the two-sided inequality 
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in Tabled 



(fe-l)r Ct Jtc 



4 /PN-^ND<S+i /p 



/ p'^^^p' S+i / p^^^nnP^' S. (38) 

For the case where homogeneous boundary conditions are prescribed on r^, i.e. 
Pn = 0, the above equations are simplified further. In such a case the total 
potential energy will coincide with the elastic strain energy and Eq. (j35p takes 
the form, 

n(t,Mc)= J Uc[j^^ccUc+^CoUt^ 

1 

In this case, the time derivative, given by Eq. p6p . is written as, 

^(i,Uc)= / Uc('^cdWd) S+i / Uo-^ddWdS. (40) 

Ot Jyc ^ JTo 

Furthermore, in this case (p^ = 0), the lower and upper energy estimates in 
the two-sided inequality are further simplified as 



Fc 

Wd.'^ddWd S, (39) 



an 



+ ^ / «-W^l)(.^OD<+^DD<-')S, (41) 
^ JFn 



and 



/ (t,u^-i)=/ u^-i(.^cD^i^-.^ce<-i)S 

J{k-i)TCt Jrc 

+ 1/" «-<-!) (.^odM^+^dd<"')S, (42) 

^ ./Fr) 



respectively. 
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The above formulation assumes the solution of each of the three sub-problems 
in all time steps, while the total response results by their superposition. It is 
worth mentioning that for proportional external loading, which is separable with 
respect to spatial coordinates and time, that means 

Uo(i, x) = ito(a;) for xeFo, (43a) 
p^{t,x) = ip{t) Pf^{x) for xgFn, (43b) 

the P° and problems need to be solved just once for the first time step, while 
for the subsequent steps the solutions will be generated by using functions (j) 
and "0 from Eq. (|43p as appropriate multipliers. 

Summarizing, expressions where the minimizer of the total potential energy 
is the displacement field defined at the adhesive contact boundary part Tc have 
been established in this section. Moreover, appropriate formulas for computing 
the lower and upper energy estimates shown in Table [2] have been given. 

4.4- Interface elements 

The interconnection of the subdomains as well as the consideration of Sig- 
norini kinematical conditions is attained by intermediate elements referred to as 
interface elements. A local reference system is associated to each interface ele- 
ment defining a normal and a tangential component of relative displacements. In 
the case of contact problems of two dcformable bodies where only small changes 
in the geometry arc assumed and conforming meshes of clastic domains are con- 
sidered along the interface, it is possible to incorporate the contact constraints 
on a purely nodal basis. For a general case of nodes arbitrarily distributed along 
the possible contact interface between two bodies, which can occur e.g. when 
automatic meshing is used for two different bodies, further considerations must 
be taken into account about the definition of Signorini contact conditions, this 
case not being considered here. The mechanical properties of springs distributed 
continuously at the interface, are given by their normal and tangential stiffnesses 
Kn and Kt , respectively, and additionally in the tangent direction also by the so- 
called plastic modulus and the factor of influence of damage kq. The shape 
functions used to approximate the distribution of variables at interface elements 
are linear and continuous for the displacements, while for the inelastic variables, 
and TT, might be constant or alternatively, continuous or discontinuous linear. 
In addition to the continuous distribution of springs, the interfaces and interface 
elements may be equipped by a "dissipative mechanism" whose properties are 
the mode-I fracture energy Gj^ and the critical stress ut, yield used in (jS]). 



5. Numerical examples 

The above introduced formulation has been implemented in a two-dimensional 
BEM code 34 1 using continuous piecewise linear boundary elements [1^, and 
also supplied with all the necessary modules for the EC-BEM, where the acronym 
EC-BEM refers to the Energetic approach for the solution of adhesive Contact 
problems by BEM. The geometry of the problem solved is shown in Figure ID 
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Figure 4; Problem geometry and boundary conditions. 

With reference to FigureO only one subdomain (i.e. = 1) is used to model in 
a simple way an experimental test motivated by the pull-push shear test used 
in engineering practice [loj . Thus, the debonding occurs between the domain 
and the rigid foundation interface Tc ■ 

The length and height of the rectangular domain fJ, respectively, are L = 250 
mm and h = 12.5 mm. The length of the initially glued part Fc placed at 
the bottom side of is Lc = 0.9L = 225 mm. The isotropic elastic mate- 
rial of the bulk is aluminium with the Young modulus E = 70 GPa and the 
Poisson ratio v ~ 0.35. Elastic plain strain state is considered. The mate- 
rial adhesive layer is epoxy resin, with elastic properties Ea = 2.4 GPa and 
Va = 0.33. Assuming the thickness of the adhesive layer ha = 0.2 mm, and 
following [4l| . the corresponding stiffness parameters are represented by the 
normal stiffness — ^ (i^i^)(i-2i/ ) GPa/mm and the tangential stiff- 

ness Kt = '''2(iZ^'')'^ The parameters for the dissipation mechanisms 

are the mode-I fracture energy Gj^ = 0.01 J/mm^ as well as the yield slip 
stress (Tt, yield =168 MPa. Then, cr„,crit = ^2KnG;^ =600MPa and crt,crit = 
y^2KtGj^ =300MPa. Finally, the hardening slope for plastic slip is Kh = Kt/9. 

5.1. Numerical experimentation 

A few sample problem cases are solved in order to illustrate the capabili- 
ties of the numerical procedure developed. In Section [5. 1.1 1 we experiment with 
a non-monotonic Dirichlet loading, for a variety of combination of dissipation 
properties. In the next Section 15.1.21 we present results for a monotonic Neu- 
mann loading on a modified geometry of the problem in order to avoid the 
absence of a Dirichlet boundary part especially after the total delamination of 
the interface. Both examples allow us to illustrate the behaviour of the numer- 
ical solution of the present energetic formulation for delamination problems, 
and do not aim to analyze the problem solutions in a thorough manner. In all 
the numerical computations, linear continuous elements have been used for the 
interface displacement and plastic slip variables, also referred to as kinematical 
variables (w, tt), while constant discontinuous elements have been assumed for 
the damage variable C- 
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Figure 5: The horizontal resultant force versus the horizontal displacement on Fd; (a) pure 
Signorini contact and (b) interface plasticity included. 



5.1.1. Non-monotonic loading 

A hard-device loading is assumed by prescribing horizontal and vertical dis- 
placements, ui{t,x) ~ sin(7t)wi(a;), where t S [0,1] while wi^l mm, and 
U2 = 0.6mi, respectively, at the right-hand side of the rectangle F, defining the 
Dirichlet boundary To. In accordance with Eg. (|43ap . (f){t)= sm{7 1) . All the 
other boundary parts are considered to be traction free, defining the Neumann 
boundary Fp,, except for the contact surface Tq. The boundary F is discretized 
by 64 elements using a uniform boundary element mesh along each side, 27 el- 
ements being used for Fc. Four combinations of properties of the dissipative 
mechanism of the adhesive are considered: 

(a) Absence of any dissipation, leading to a pure elastic Signorini contact prob- 
lem, 

(b) Interface plasticity is considered, the damage variables C, being excluded 
from the minimization procedure, 

(c) Interface damage is considered, the plastic slip variables tt being excluded 
from the minimization procedure, 

(d) Both interface damage and plasticity are considered. 

Cases (a) and (b) are mainly included for the comparison purposes and also in 
order to analyse an inelastic response due to interface plasticity. The horizontal 
resultant force with respect to the displacement on F^ is plotted in Fig. [S] For 
the case of an inelastic response due to a non-monotonic loading, a hysteresis 
cycle appears as expected. Furthermore, for these two cases also the shear 
stresses with respect to the relative tangential displacements at an interface 
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Figure 6: Stress-relative displacement behaviour computed at the interface point xi = 208.33 
mm of Fc : (a) pure Signorini contact and (b) interface plasticity included. 



point Xi = 208.33 mm are plotted in Fig. [SI where a typical hysteretic behaviour 
for the kinematic type hardening plasticity is successfully computed in case (b) . 

More complicated behaviour is obtained for cases (c) and (d) where inter- 
face damage is included. This may be observed in Figs. [7] and |8] where the 
horizontal and vertical resultant forces with the respective displacements at Fd 
are depicted. For these cases upon the first uploading a damage initially ap- 
pears, for case (c) by breaking one element that corresponds to a crack opening 
of Lcrack = 0.037ic = 8.33 mm, while for case (d) in a following time step, by 
breaking six elements simultaneously which corresponds to a crack opening of 
^crack = 0.22Lc = 50.0 mm. For the same time step, that crack initiates in case 
(d) with a crack opening of icrack — 50.0 mm, for case (c) after some progressive 
damage propagation, a crack of the same length exists. This behaviour is the 
expected one, since because of plasticity appearance in case (d), damage delayed 
to appear in comparison with case c) where it is assumed that an energy may 
be released only due to damage. Then, after change of the direction of load- 
ing no further damage appears, while plasticity still evolves on the remaining 
glued part of Fc upon the respective uploading periods. This behavior can be 
better understood from Fig. |9l where the evolution of the accumulated dissipa- 
tion with respect to the time t is shown. In fact, t is a kind of pseudo-time or 
process time which can arbitrarily be re-scaled since the considered system is 
rate- independent. Finally, the evolution of stored energies in the adhesive layer 
due to opening and shear are shown in Fig. 1101 
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Figure 7: The horizontal resultant force versus the horizontal displacement on Fd : (c) interface 
damage and (d) both interface damage and plasticity. 




Figure 8: The vertical resultant force versus the vertical displacement on Fd: (c) interface 
damage and (d) both interface damage and plasticity. 
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Figure 9: Evolution of the dissipated energies: (c) interface damage and (d) both interface 
damage and plasticity. 
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Figure 10: Evolution of the stored energies in the adhesive for case (d) considering both 
interface damage and plasticity. 
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Figure 11: Modified geometry and boundary conditions used for the problem with prescribed 
non-zero tractions. 
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Figure 12; The horizontal resultant force on Fiin versus the horizontal displacement of the 
lower right corner of H. 



5.1.2. Traction instead of displacement loading 

From an engineering point of view, we are highly interested in the case of 
external loading given by nonvanishing Neumann boundary conditions. This 
is because in numerous experiments or real applications, loading is described 
through external forces, moments or tractions. For these reasons a modified 
problem configuration shown in Fig. [TT] is studied in this section. The length of 
Fc is Lc = 200 mm, while the homogeneous Neumann boundary parts F^/i on 
the left and right hand side of Fc have lengths equal to 0.125Lc. The length and 
height of O as well as the material properties are the same as in the previous 
example. The left vertical side of D, is fixed, defining Fo, while uniform normal 
tractions applied on the right vertical side of 51, defining Fnti, arc increasing in 
time, i.e. a;) = 'ip{t)po and P2 = therein, with pq > being a constant 

and ip{t) > an increasing function, see Eq. (j43b|. Both interface damage and 
plasticity are considered. 

The evolution of the horizontal resultant force on Fn„ is plotted in Fig. [T2l 
The computational analysis, including BTA together with the two-sided energy 
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Figure 13: Shear stress-relative tangential displacement relation as computed at the rightmost 
point of the interface Fc. For this case, only plasticity is considered, damage being excluded. 



inequality checking, stops at a point, marked in Fig. [T^ as point (L), where after 
some plasticity development the first and at the same time the total damage of 
the adhesive layer occurs. A dashed line in the same plot represents the tangent 
line to the initial purely elastic part of the resultant force-displacement curve. 
Both lines separate due to the appearance of some plastic slip at Tq- To capture 
a progressive damage propagation along Fc would require decreasing the applied 
load after the peak load is achieved. 

Finally, in order to analyse the pointwise behaviour of the numerical solu- 
tion at Fc, in particular regarding the evolution of plastic slip, the same problem 
is solved again but including interface plasticity only, i.e. no damage at Fq is 
possible. Numerically computed shear stress versus the relative tangential dis- 
placement at the rightmost node of Fc is shown in Fig. 1131 together with the 
expected tangential stress-relative displacement law of Fig. [5J An overshoot- 
ing^^ phenomenon takes place when plasticity occurs, a similar behaviour may 
also be observed in Fig. |6l This phenomenon is essentially associated to the time 
and spatial discretization of the problem, in particular possibly due to some os- 
cillations of the traction solution near the crack tip, and therefore can gradually 
be eliminated by a spatial discretization refinement as observed in Fig. 1131 

5.2. Practical application 

The problem configuration shown in Fig. 2] is considered again. A mono- 
tonic hard-device loading is assumed by prescribing horizontal and vertical 
displacements, respectively, as ui{t,x) = twi{x) with wi{x) = 0.6mm and 
U2it,x) = 0.6ui(t, x) at the right-hand side of the rectangle F, defining the 
Dirichlet boundary Fq. All the other boundary parts are considered to be trac- 
tion free, defining the Neumann boundary F^, except for the contact surface Fc. 
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The problem evolution is represented as function of the fictitious time t. 

An advantage of the present method considering a continuous distribution of 
springs along the interface in comparison with the classical fracture mechanics, 
which assumes a perfect interface except for some cracked zones, is that no 
special mesh refinement is needed near the crack tip and a uniform mesh can be 
employed along Fc, similarly as in the Cohesive Zone Models, see 18| and further 
references therein. The tractions along the adhesive layer of the present type 
are bounded although traction concentrations can be expected at the end-points 



of the adhesive layer, which may correspond to crack tips, cf. |23|,|40|. Actually, 
in the present model, these tractions are limited by the critical values of normal 
and tangential tractions cr„^crit and trt^crit, respectively. If the adhesive layer 
at Fc in the present problem in Fig. |4] would be replaced by perfect bonding 
conditions, stress singularities would appear at both extremes of the bonded 
part Fc. The stress singularity at the right extreme of Fc, corresponding to the 
classical oscillatory singularity of the open model of an interface crack between 
an elastic and an infinitely rigid solid, cf. (4^ . would be more severe than that 
at the left extreme. In such strongly refined mesh or special singular 

elements would be needed for a proper problem discretization of the crack tip 
neighbourhood. It should be mentioned here that intuitive refinement without 
having rigorous local error indicators is sometimes dangerous and may destroy 
convergence which is standardly guaranteed on uniformly refined meshes only. 

Nevertheless, in the present model, the fact that a local mesh refinement 
at the crack tip is not needed makes easy the modeling of damage progression 
with the crack tip moving along the interface. In order to check this statement, 
we have solved the problem in Fig. |4] by using a series of uniform boundary 
element meshes, the three finest meshes having 126 (60 and 3 elements along 
each horizontal and vertical side, respectively) , 252 and 504 elements on F (that 
corresponds to 54, 108 and 216 elements on Fc). The traction solutions along Fc 
for these three finest meshes shown in Fig.[T4]correspond to horizontal prescribed 
displacement Ui = 0.28mm, when no damage appears although some interface 
plasticity has evolved. A strong traction concentration at the right extreme of 
Fc can be observed in these plots which, however, indicate that even for the 
coarsest mesh (54 elements on Fc) the solution obtained is sufiiciently accurate 
for the purpose of the present study. An additional checking is shown in Fig. I15i 
where the percentage differences of the computed strain energy in the bulk, the 
computed total energy (that is the sum of the stored energy and the dissipated 
energy at time t, S'{t,u{t), z{t)) + Diss-R,(z; [0, t], see ([2])), and the maximum 
absolute value of normal tractions at Fc computed by a coarse mesh and the 
finest mesh (216 elements on Fc) are plotted. These plots confirm that the 
percentage difference of the strain energy and of the maximum normal traction 
for the mesh with 54 elements on Fc is sufficiently small, in particular it is about 
1%. Therefore, this mesh is used in the following complete numerical study of 
the present problem. 

Fig. [16] shows the evolution of different energies computed, in particular, 
the energy stored in the elastic bulk, in the adhesive layer and the dissipated 
energy. Also the total energy, which is actually minimized in the time stepping 
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Figure 14: (a) Normal and (b) tangent tractions along the adhesive zone Fc, for the three 
finest uniform boundary element meshes, n the number of boundary elements at Fc . 
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Figure 15: Percentage difference of (a) tlie total energy, (b) the elastic strain energy and (c) 
the maximum normal tractions at the adhesive zone Fc, for uniform, along the horizontal and 
vertical edges, boundary element meshes. 
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Figure 17: Deformed shape of the elastic domain for two time steps, just before and after 
damage initiation. A scale factor of 150 has been used. 



procedure, together with the lower and upper estimates of energy are shown. As 
it can be seen in Fig. 1161 the global minimization procedure defines the end of 
the delamination process, where the remaining undamaged part of the adhesive 
layer is dcbonded, as the point where the sum of the stored energy in the bulk 
and the adhesive layer is basically equal to the energy needed to delaminate the 
undamaged part of the adhesive layer, which is given by the dissipated energy 
in the last time step. In Fig. [TTl the deformed shape of the bulk, is plotted for 
two time steps, just before and after the first damage, respectively. 

Figs. [18] and [19] present the components of normal and tangential traction 
vector along the adhesive zone Fc. A very good agreement exists between the 
computed tractions for each subdomain by BEM and those computed in the 
adhesive layer, although the equilibrium has not been imposed directly but it 
results as a consequence of the energy minimization. Progressive extension of 
the traction free portion of the original Fc because of the damage propagation 
(C ~ 0) can be observed in Figs. [T8b) and (TOb) . It should be mentioned that 
the portion of Fc which is totally damaged is still kept as a part of the mini- 
mization procedure, where nodal displacement values participate as unknowns 
in the minimization procedure and their values are used in the BEM solution 
of the pertinent BVP. For this reason, in fact an approximation of the devel- 
oped traction-free zone is computed by BEM for each subdomain. Obviously 
other algorithms might be used where after the total damage of a portion of the 
adhesive layer a switch in the type of the boundary condition (e.g. from Dirich- 
let to vanishing Neumann boundary condition) is taken into account along this 
boundary portion in the BEM computation for each subdomain. Nevertheless, 
we have been interested in the results obtained by the present simple procedure. 
Normal compressive tractions computed by BEM can be observed in Fig. 1181 in 
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zones where vanishing normal tractions are obtained in the interface elements. 
This is due to the fact that the rigid obstacle undertakes these compressive 
tractions. 

Finally, in Fig. [50] the resultant forces acting at versus the pseudo-time 
t are shown. These two plots have some similarities in the behaviour that 
may come up through the characteristic points (A)-(E). Up to point (A) the 
linear elastic behaviour manifests in both the solid and adhesive, at this point a 
plastic slip appears. Then, at point (B) the first damage appears in the first 11 
boundary elements that are situated on the right-hand side of the adhesive layer. 
This new crack length results in a "jump down" of the resultant forces up to 
point (C). Then, up to point (D) the damaged zone is progressively extended and 
finally after point (D) the remaining adhesive zone is damaged instantaneously. 
The problem evolution ends up at point (E), where a rigid body motion of the 
elastic body takes place. The increment of the crack length from point (B) 
to (C) equals 45.83 mm. In the same figures also the linear elastic responses, 
obtained for the same configuration taking into account only Signorini contact 
without any interface damage and plasticity, are plotted in order to facilitate 
the observation of the initiation of plasticity and/or damage. 

6. Conclusions 

A boundary element implementation of a computational procedure based on 
an energetic-solution framework for the delamination problems has been pre- 
sented. A specific model for the adhesive interfaces, which distinguishes the 
amount of energy dissipated in opening Mode I and shear Mode II has been 
adopted. This model involves two inelastic internal variables on delaminating 
surfaces, namely the damage variable C and the plastic slip tt. Some details 
regarding the formulation of the collocation BEM as well as the optimization 
procedures necessary for solving the global minimization problem, inherent in 
the formulation, have been discussed. A few numerical tests have been pre- 
sented in order to analyse the behaviour of the present delamination model and 
performance of the algorithms implemented in a collocational BEM code. 
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Figure 18: Normal tractions along the adhesive zone Fc, just (a) before and (b) after the first 
crack opening, computed by BEM as well as in the interface elements. 
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Figure 20: Evolution of (a) the horizontal and (b) vertical resultant force on 
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